using System;

namespace ICP_MS.Plugin.HPLC;

public class Matrix
{
	private int numColumns = 0;

	private int numRows = 0;

	private double eps = 0.0;

	private double[] elements = null;

	public int Columns => numColumns;

	public int Rows => numRows;

	public double this[int row, int col]
	{
		get
		{
			return elements[col + row * numColumns];
		}
		set
		{
			elements[col + row * numColumns] = value;
		}
	}

	public double Eps
	{
		get
		{
			return eps;
		}
		set
		{
			eps = value;
		}
	}

	public Matrix()
	{
		numColumns = 1;
		numRows = 1;
		Init(numRows, numColumns);
	}

	public Matrix(int nRows, int nCols)
	{
		numRows = nRows;
		numColumns = nCols;
		Init(numRows, numColumns);
	}

	public Matrix(double[,] value)
	{
		numRows = value.GetLength(0);
		numColumns = value.GetLength(1);
		double[] data = new double[numRows * numColumns];
		int k = 0;
		for (int i = 0; i < numRows; i++)
		{
			for (int j = 0; j < numColumns; j++)
			{
				data[k++] = value[i, j];
			}
		}
		Init(numRows, numColumns);
		SetData(data);
	}

	public Matrix(int nRows, int nCols, double[] value)
	{
		numRows = nRows;
		numColumns = nCols;
		Init(numRows, numColumns);
		SetData(value);
	}

	public Matrix(int nSize)
	{
		numRows = nSize;
		numColumns = nSize;
		Init(nSize, nSize);
	}

	public Matrix(int nSize, double[] value)
	{
		numRows = nSize;
		numColumns = nSize;
		Init(nSize, nSize);
		SetData(value);
	}

	public Matrix(Matrix other)
	{
		numColumns = other.GetNumColumns();
		numRows = other.GetNumRows();
		Init(numRows, numColumns);
		SetData(other.elements);
	}

	public bool Init(int nRows, int nCols)
	{
		numRows = nRows;
		numColumns = nCols;
		int nSize = nCols * nRows;
		if (nSize < 0)
		{
			return false;
		}
		elements = new double[nSize];
		return true;
	}

	public void SetEps(double newEps)
	{
		eps = newEps;
	}

	public double GetEps()
	{
		return eps;
	}

	public static Matrix operator +(Matrix m1, Matrix m2)
	{
		return m1.Add(m2);
	}

	public static Matrix operator -(Matrix m1, Matrix m2)
	{
		return m1.Subtract(m2);
	}

	public static Matrix operator *(Matrix m1, Matrix m2)
	{
		return m1.Multiply(m2);
	}

	public static implicit operator double[](Matrix m)
	{
		return m.elements;
	}

	public bool MakeUnitMatrix(int nSize)
	{
		if (!Init(nSize, nSize))
		{
			return false;
		}
		for (int i = 0; i < nSize; i++)
		{
			for (int j = 0; j < nSize; j++)
			{
				if (i == j)
				{
					SetElement(i, j, 1.0);
				}
			}
		}
		return true;
	}

	public override string ToString()
	{
		return ToString(",", bLineBreak: true);
	}

	public string ToString(string sDelim, bool bLineBreak)
	{
		string s = "";
		for (int i = 0; i < numRows; i++)
		{
			for (int j = 0; j < numColumns; j++)
			{
				string ss = GetElement(i, j).ToString("F");
				s += ss;
				if (bLineBreak)
				{
					if (j != numColumns - 1)
					{
						s += sDelim;
					}
				}
				else if (i != numRows - 1 || j != numColumns - 1)
				{
					s += sDelim;
				}
			}
			if (bLineBreak && i != numRows - 1)
			{
				s += "\r\n";
			}
		}
		return s;
	}

	public string ToStringRow(int nRow, string sDelim)
	{
		string s = "";
		if (nRow >= numRows)
		{
			return s;
		}
		for (int i = 0; i < numColumns; i++)
		{
			string ss = GetElement(nRow, i).ToString("F");
			s += ss;
			if (i != numColumns - 1)
			{
				s += sDelim;
			}
		}
		return s;
	}

	public string ToStringCol(int nCol, string sDelim)
	{
		string s = "";
		if (nCol >= numColumns)
		{
			return s;
		}
		for (int i = 0; i < numRows; i++)
		{
			string ss = GetElement(i, nCol).ToString("F");
			s += ss;
			if (i != numRows - 1)
			{
				s += sDelim;
			}
		}
		return s;
	}

	public void SetData(double[] value)
	{
		elements = (double[])value.Clone();
	}

	public bool SetElement(int nRow, int nCol, double value)
	{
		if (nCol < 0 || nCol >= numColumns || nRow < 0 || nRow >= numRows)
		{
			return false;
		}
		elements[nCol + nRow * numColumns] = value;
		return true;
	}

	public double GetElement(int nRow, int nCol)
	{
		return elements[nCol + nRow * numColumns];
	}

	public int GetNumColumns()
	{
		return numColumns;
	}

	public int GetNumRows()
	{
		return numRows;
	}

	public double[] GetData()
	{
		return elements;
	}

	public int GetRowVector(int nRow, double[] pVector)
	{
		for (int i = 0; i < numColumns; i++)
		{
			pVector[i] = GetElement(nRow, i);
		}
		return numColumns;
	}

	public int GetColVector(int nCol, double[] pVector)
	{
		for (int i = 0; i < numRows; i++)
		{
			pVector[i] = GetElement(i, nCol);
		}
		return numRows;
	}

	public Matrix SetValue(Matrix other)
	{
		if (other != this)
		{
			Init(other.GetNumRows(), other.GetNumColumns());
			SetData(other.elements);
		}
		return this;
	}

	public override bool Equals(object other)
	{
		if (!(other is Matrix matrix))
		{
			return false;
		}
		if (numColumns != matrix.GetNumColumns() || numRows != matrix.GetNumRows())
		{
			return false;
		}
		for (int i = 0; i < numRows; i++)
		{
			for (int j = 0; j < numColumns; j++)
			{
				if (Math.Abs(GetElement(i, j) - matrix.GetElement(i, j)) > eps)
				{
					return false;
				}
			}
		}
		return true;
	}

	public override int GetHashCode()
	{
		double sum = 0.0;
		for (int i = 0; i < numRows; i++)
		{
			for (int j = 0; j < numColumns; j++)
			{
				sum += Math.Abs(GetElement(i, j));
			}
		}
		return (int)Math.Sqrt(sum);
	}

	public Matrix Add(Matrix other)
	{
		if (numColumns != other.GetNumColumns() || numRows != other.GetNumRows())
		{
			throw new Exception("矩阵的行/列数不匹配。");
		}
		Matrix result = new Matrix(this);
		for (int i = 0; i < numRows; i++)
		{
			for (int j = 0; j < numColumns; j++)
			{
				result.SetElement(i, j, result.GetElement(i, j) + other.GetElement(i, j));
			}
		}
		return result;
	}

	public Matrix Subtract(Matrix other)
	{
		if (numColumns != other.GetNumColumns() || numRows != other.GetNumRows())
		{
			throw new Exception("矩阵的行/列数不匹配。");
		}
		Matrix result = new Matrix(this);
		for (int i = 0; i < numRows; i++)
		{
			for (int j = 0; j < numColumns; j++)
			{
				result.SetElement(i, j, result.GetElement(i, j) - other.GetElement(i, j));
			}
		}
		return result;
	}

	public Matrix Multiply(double value)
	{
		Matrix result = new Matrix(this);
		for (int i = 0; i < numRows; i++)
		{
			for (int j = 0; j < numColumns; j++)
			{
				result.SetElement(i, j, result.GetElement(i, j) * value);
			}
		}
		return result;
	}

	public Matrix Multiply(Matrix other)
	{
		if (numColumns != other.GetNumRows())
		{
			throw new Exception("矩阵的行/列数不匹配。");
		}
		Matrix result = new Matrix(numRows, other.GetNumColumns());
		for (int i = 0; i < result.GetNumRows(); i++)
		{
			for (int j = 0; j < other.GetNumColumns(); j++)
			{
				double value = 0.0;
				for (int k = 0; k < numColumns; k++)
				{
					value += GetElement(i, k) * other.GetElement(k, j);
				}
				result.SetElement(i, j, value);
			}
		}
		return result;
	}

	public bool Multiply(Matrix AR, Matrix AI, Matrix BR, Matrix BI, Matrix CR, Matrix CI)
	{
		if (AR.GetNumColumns() != AI.GetNumColumns() || AR.GetNumRows() != AI.GetNumRows() || BR.GetNumColumns() != BI.GetNumColumns() || BR.GetNumRows() != BI.GetNumRows() || AR.GetNumColumns() != BR.GetNumRows())
		{
			return false;
		}
		Matrix mtxCR = new Matrix(AR.GetNumRows(), BR.GetNumColumns());
		Matrix mtxCI = new Matrix(AR.GetNumRows(), BR.GetNumColumns());
		for (int i = 0; i < AR.GetNumRows(); i++)
		{
			for (int j = 0; j < BR.GetNumColumns(); j++)
			{
				double vr = 0.0;
				double vi = 0.0;
				for (int k = 0; k < AR.GetNumColumns(); k++)
				{
					double p = AR.GetElement(i, k) * BR.GetElement(k, j);
					double q = AI.GetElement(i, k) * BI.GetElement(k, j);
					double s = (AR.GetElement(i, k) + AI.GetElement(i, k)) * (BR.GetElement(k, j) + BI.GetElement(k, j));
					vr += p - q;
					vi += s - p - q;
				}
				mtxCR.SetElement(i, j, vr);
				mtxCI.SetElement(i, j, vi);
			}
		}
		CR = mtxCR;
		CI = mtxCI;
		return true;
	}

	public Matrix Transpose()
	{
		Matrix Trans = new Matrix(numColumns, numRows);
		for (int i = 0; i < numRows; i++)
		{
			for (int j = 0; j < numColumns; j++)
			{
				Trans.SetElement(j, i, GetElement(i, j));
			}
		}
		return Trans;
	}

	public bool InvertGaussJordan()
	{
		double d = 0.0;
		double p = 0.0;
		int[] pnRow = new int[numColumns];
		int[] pnCol = new int[numColumns];
		for (int k = 0; k <= numColumns - 1; k++)
		{
			d = 0.0;
			int l;
			for (int i = k; i <= numColumns - 1; i++)
			{
				for (int j = k; j <= numColumns - 1; j++)
				{
					l = i * numColumns + j;
					p = Math.Abs(elements[l]);
					if (p > d)
					{
						d = p;
						pnRow[k] = i;
						pnCol[k] = j;
					}
				}
			}
			if (d == 0.0)
			{
				return false;
			}
			if (pnRow[k] != k)
			{
				for (int j = 0; j <= numColumns - 1; j++)
				{
					int u = k * numColumns + j;
					int v = pnRow[k] * numColumns + j;
					p = elements[u];
					elements[u] = elements[v];
					elements[v] = p;
				}
			}
			if (pnCol[k] != k)
			{
				for (int i = 0; i <= numColumns - 1; i++)
				{
					int u = i * numColumns + k;
					int v = i * numColumns + pnCol[k];
					p = elements[u];
					elements[u] = elements[v];
					elements[v] = p;
				}
			}
			l = k * numColumns + k;
			elements[l] = 1.0 / elements[l];
			for (int j = 0; j <= numColumns - 1; j++)
			{
				if (j != k)
				{
					int u = k * numColumns + j;
					elements[u] *= elements[l];
				}
			}
			for (int i = 0; i <= numColumns - 1; i++)
			{
				if (i == k)
				{
					continue;
				}
				for (int j = 0; j <= numColumns - 1; j++)
				{
					if (j != k)
					{
						int u = i * numColumns + j;
						elements[u] -= elements[i * numColumns + k] * elements[k * numColumns + j];
					}
				}
			}
			for (int i = 0; i <= numColumns - 1; i++)
			{
				if (i != k)
				{
					int u = i * numColumns + k;
					elements[u] = (0.0 - elements[u]) * elements[l];
				}
			}
		}
		for (int k = numColumns - 1; k >= 0; k--)
		{
			if (pnCol[k] != k)
			{
				for (int j = 0; j <= numColumns - 1; j++)
				{
					int u = k * numColumns + j;
					int v = pnCol[k] * numColumns + j;
					p = elements[u];
					elements[u] = elements[v];
					elements[v] = p;
				}
			}
			if (pnRow[k] != k)
			{
				for (int i = 0; i <= numColumns - 1; i++)
				{
					int u = i * numColumns + k;
					int v = i * numColumns + pnRow[k];
					p = elements[u];
					elements[u] = elements[v];
					elements[v] = p;
				}
			}
		}
		return true;
	}

	public bool InvertGaussJordan(Matrix mtxImag)
	{
		int[] pnRow = new int[numColumns];
		int[] pnCol = new int[numColumns];
		for (int k = 0; k <= numColumns - 1; k++)
		{
			double d = 0.0;
			for (int i = k; i <= numColumns - 1; i++)
			{
				for (int j = k; j <= numColumns - 1; j++)
				{
					int u = i * numColumns + j;
					double p = elements[u] * elements[u] + mtxImag.elements[u] * mtxImag.elements[u];
					if (p > d)
					{
						d = p;
						pnRow[k] = i;
						pnCol[k] = j;
					}
				}
			}
			if (d == 0.0)
			{
				return false;
			}
			if (pnRow[k] != k)
			{
				for (int j = 0; j <= numColumns - 1; j++)
				{
					int u = k * numColumns + j;
					int v = pnRow[k] * numColumns + j;
					double t = elements[u];
					elements[u] = elements[v];
					elements[v] = t;
					t = mtxImag.elements[u];
					mtxImag.elements[u] = mtxImag.elements[v];
					mtxImag.elements[v] = t;
				}
			}
			if (pnCol[k] != k)
			{
				for (int i = 0; i <= numColumns - 1; i++)
				{
					int u = i * numColumns + k;
					int v = i * numColumns + pnCol[k];
					double t = elements[u];
					elements[u] = elements[v];
					elements[v] = t;
					t = mtxImag.elements[u];
					mtxImag.elements[u] = mtxImag.elements[v];
					mtxImag.elements[v] = t;
				}
			}
			int l = k * numColumns + k;
			elements[l] /= d;
			mtxImag.elements[l] = (0.0 - mtxImag.elements[l]) / d;
			for (int j = 0; j <= numColumns - 1; j++)
			{
				if (j != k)
				{
					int u = k * numColumns + j;
					double p = elements[u] * elements[l];
					double q = mtxImag.elements[u] * mtxImag.elements[l];
					double s = (elements[u] + mtxImag.elements[u]) * (elements[l] + mtxImag.elements[l]);
					elements[u] = p - q;
					mtxImag.elements[u] = s - p - q;
				}
			}
			for (int i = 0; i <= numColumns - 1; i++)
			{
				if (i == k)
				{
					continue;
				}
				int v = i * numColumns + k;
				for (int j = 0; j <= numColumns - 1; j++)
				{
					if (j != k)
					{
						int u = k * numColumns + j;
						int w = i * numColumns + j;
						double p = elements[u] * elements[v];
						double q = mtxImag.elements[u] * mtxImag.elements[v];
						double s = (elements[u] + mtxImag.elements[u]) * (elements[v] + mtxImag.elements[v]);
						double t = p - q;
						double b = s - p - q;
						elements[w] -= t;
						mtxImag.elements[w] = mtxImag.elements[w] - b;
					}
				}
			}
			for (int i = 0; i <= numColumns - 1; i++)
			{
				if (i != k)
				{
					int u = i * numColumns + k;
					double p = elements[u] * elements[l];
					double q = mtxImag.elements[u] * mtxImag.elements[l];
					double s = (elements[u] + mtxImag.elements[u]) * (elements[l] + mtxImag.elements[l]);
					elements[u] = q - p;
					mtxImag.elements[u] = p + q - s;
				}
			}
		}
		for (int k = numColumns - 1; k >= 0; k--)
		{
			if (pnCol[k] != k)
			{
				for (int j = 0; j <= numColumns - 1; j++)
				{
					int u = k * numColumns + j;
					int v = pnCol[k] * numColumns + j;
					double t = elements[u];
					elements[u] = elements[v];
					elements[v] = t;
					t = mtxImag.elements[u];
					mtxImag.elements[u] = mtxImag.elements[v];
					mtxImag.elements[v] = t;
				}
			}
			if (pnRow[k] != k)
			{
				for (int i = 0; i <= numColumns - 1; i++)
				{
					int u = i * numColumns + k;
					int v = i * numColumns + pnRow[k];
					double t = elements[u];
					elements[u] = elements[v];
					elements[v] = t;
					t = mtxImag.elements[u];
					mtxImag.elements[u] = mtxImag.elements[v];
					mtxImag.elements[v] = t;
				}
			}
		}
		return true;
	}

	public bool InvertSsgj()
	{
		double[] pTmp = new double[numColumns];
		for (int k = 0; k <= numColumns - 1; k++)
		{
			double w = elements[0];
			if (w == 0.0)
			{
				return false;
			}
			int l = numColumns - k - 1;
			for (int i = 1; i <= numColumns - 1; i++)
			{
				double g = elements[i * numColumns];
				pTmp[i] = g / w;
				if (i <= l)
				{
					pTmp[i] = 0.0 - pTmp[i];
				}
				for (int j = 1; j <= i; j++)
				{
					elements[(i - 1) * numColumns + j - 1] = elements[i * numColumns + j] + g * pTmp[j];
				}
			}
			elements[numColumns * numColumns - 1] = 1.0 / w;
			for (int i = 1; i <= numColumns - 1; i++)
			{
				elements[(numColumns - 1) * numColumns + i - 1] = pTmp[i];
			}
		}
		for (int i = 0; i <= numColumns - 2; i++)
		{
			for (int j = i + 1; j <= numColumns - 1; j++)
			{
				elements[i * numColumns + j] = elements[j * numColumns + i];
			}
		}
		return true;
	}

	public bool InvertTrench()
	{
		double[] t = new double[numColumns];
		double[] tt = new double[numColumns];
		for (int i = 0; i < numColumns; i++)
		{
			t[i] = GetElement(0, i);
			tt[i] = GetElement(i, 0);
		}
		double[] c = new double[numColumns];
		double[] r = new double[numColumns];
		double[] p = new double[numColumns];
		if (t[0] == 0.0)
		{
			return false;
		}
		double a = t[0];
		c[0] = tt[1] / t[0];
		r[0] = t[1] / t[0];
		for (int k = 0; k <= numColumns - 3; k++)
		{
			double s = 0.0;
			for (int j = 1; j <= k + 1; j++)
			{
				s += c[k + 1 - j] * tt[j];
			}
			s = (s - tt[k + 2]) / a;
			for (int i = 0; i <= k; i++)
			{
				p[i] = c[i] + s * r[k - i];
			}
			c[k + 1] = 0.0 - s;
			s = 0.0;
			for (int j = 1; j <= k + 1; j++)
			{
				s += r[k + 1 - j] * t[j];
			}
			s = (s - t[k + 2]) / a;
			for (int i = 0; i <= k; i++)
			{
				r[i] += s * c[k - i];
				c[k - i] = p[k - i];
			}
			r[k + 1] = 0.0 - s;
			a = 0.0;
			for (int j = 1; j <= k + 2; j++)
			{
				a += t[j] * c[j - 1];
			}
			a = t[0] - a;
			if (a == 0.0)
			{
				return false;
			}
		}
		elements[0] = 1.0 / a;
		for (int i = 0; i <= numColumns - 2; i++)
		{
			int k = i + 1;
			int j = (i + 1) * numColumns;
			elements[k] = (0.0 - r[i]) / a;
			elements[j] = (0.0 - c[i]) / a;
		}
		for (int i = 0; i <= numColumns - 2; i++)
		{
			for (int j = 0; j <= numColumns - 2; j++)
			{
				int k = (i + 1) * numColumns + j + 1;
				elements[k] = elements[i * numColumns + j] - c[i] * elements[j + 1];
				elements[k] += c[numColumns - j - 2] * elements[numColumns - i - 1];
			}
		}
		return true;
	}

	public double ComputeDetGauss()
	{
		int nis = 0;
		int js = 0;
		double f = 1.0;
		double det = 1.0;
		for (int k = 0; k <= numColumns - 2; k++)
		{
			double q = 0.0;
			int l;
			for (int i = k; i <= numColumns - 1; i++)
			{
				for (int j = k; j <= numColumns - 1; j++)
				{
					l = i * numColumns + j;
					double d = Math.Abs(elements[l]);
					if (d > q)
					{
						q = d;
						nis = i;
						js = j;
					}
				}
			}
			if (q == 0.0)
			{
				return 0.0;
			}
			if (nis != k)
			{
				f = 0.0 - f;
				for (int j = k; j <= numColumns - 1; j++)
				{
					int u = k * numColumns + j;
					int v = nis * numColumns + j;
					double d = elements[u];
					elements[u] = elements[v];
					elements[v] = d;
				}
			}
			if (js != k)
			{
				f = 0.0 - f;
				for (int i = k; i <= numColumns - 1; i++)
				{
					int u = i * numColumns + js;
					int v = i * numColumns + k;
					double d = elements[u];
					elements[u] = elements[v];
					elements[v] = d;
				}
			}
			l = k * numColumns + k;
			det *= elements[l];
			for (int i = k + 1; i <= numColumns - 1; i++)
			{
				double d = elements[i * numColumns + k] / elements[l];
				for (int j = k + 1; j <= numColumns - 1; j++)
				{
					int u = i * numColumns + j;
					elements[u] -= d * elements[k * numColumns + j];
				}
			}
		}
		return f * det * elements[numColumns * numColumns - 1];
	}

	public int ComputeRankGauss()
	{
		int nis = 0;
		int js = 0;
		int nn = numRows;
		if (numRows >= numColumns)
		{
			nn = numColumns;
		}
		int k = 0;
		for (int l = 0; l <= nn - 1; l++)
		{
			double q = 0.0;
			int ll;
			for (int i = l; i <= numRows - 1; i++)
			{
				for (int j = l; j <= numColumns - 1; j++)
				{
					ll = i * numColumns + j;
					double d = Math.Abs(elements[ll]);
					if (d > q)
					{
						q = d;
						nis = i;
						js = j;
					}
				}
			}
			if (q == 0.0)
			{
				return k;
			}
			k++;
			if (nis != l)
			{
				for (int j = l; j <= numColumns - 1; j++)
				{
					int u = l * numColumns + j;
					int v = nis * numColumns + j;
					double d = elements[u];
					elements[u] = elements[v];
					elements[v] = d;
				}
			}
			if (js != l)
			{
				for (int i = l; i <= numRows - 1; i++)
				{
					int u = i * numColumns + js;
					int v = i * numColumns + l;
					double d = elements[u];
					elements[u] = elements[v];
					elements[v] = d;
				}
			}
			ll = l * numColumns + l;
			for (int i = l + 1; i <= numColumns - 1; i++)
			{
				double d = elements[i * numColumns + l] / elements[ll];
				for (int j = l + 1; j <= numColumns - 1; j++)
				{
					int u = i * numColumns + j;
					elements[u] -= d * elements[l * numColumns + j];
				}
			}
		}
		return k;
	}

	public bool ComputeDetCholesky(ref double realDetValue)
	{
		if (elements[0] <= 0.0)
		{
			return false;
		}
		elements[0] = Math.Sqrt(elements[0]);
		double d = elements[0];
		for (int i = 1; i <= numColumns - 1; i++)
		{
			int u = i * numColumns;
			elements[u] /= elements[0];
		}
		for (int j = 1; j <= numColumns - 1; j++)
		{
			int l = j * numColumns + j;
			for (int k = 0; k <= j - 1; k++)
			{
				int u = j * numColumns + k;
				elements[l] -= elements[u] * elements[u];
			}
			if (elements[l] <= 0.0)
			{
				return false;
			}
			elements[l] = Math.Sqrt(elements[l]);
			d *= elements[l];
			for (int i = j + 1; i <= numColumns - 1; i++)
			{
				int u = i * numColumns + j;
				for (int k = 0; k <= j - 1; k++)
				{
					elements[u] -= elements[i * numColumns + k] * elements[j * numColumns + k];
				}
				elements[u] /= elements[l];
			}
		}
		realDetValue = d * d;
		for (int i = 0; i <= numColumns - 2; i++)
		{
			for (int j = i + 1; j <= numColumns - 1; j++)
			{
				elements[i * numColumns + j] = 0.0;
			}
		}
		return true;
	}

	public bool SplitLU(Matrix mtxL, Matrix mtxU)
	{
		if (!mtxL.Init(numColumns, numColumns) || !mtxU.Init(numColumns, numColumns))
		{
			return false;
		}
		for (int k = 0; k <= numColumns - 2; k++)
		{
			int ll = k * numColumns + k;
			if (elements[ll] == 0.0)
			{
				return false;
			}
			for (int i = k + 1; i <= numColumns - 1; i++)
			{
				int w = i * numColumns + k;
				elements[w] /= elements[ll];
			}
			for (int i = k + 1; i <= numColumns - 1; i++)
			{
				int w = i * numColumns + k;
				for (int j = k + 1; j <= numColumns - 1; j++)
				{
					int v = i * numColumns + j;
					elements[v] -= elements[w] * elements[k * numColumns + j];
				}
			}
		}
		for (int i = 0; i <= numColumns - 1; i++)
		{
			int w;
			for (int j = 0; j < i; j++)
			{
				w = i * numColumns + j;
				mtxL.elements[w] = elements[w];
				mtxU.elements[w] = 0.0;
			}
			w = i * numColumns + i;
			mtxL.elements[w] = 1.0;
			mtxU.elements[w] = elements[w];
			for (int j = i + 1; j <= numColumns - 1; j++)
			{
				w = i * numColumns + j;
				mtxL.elements[w] = 0.0;
				mtxU.elements[w] = elements[w];
			}
		}
		return true;
	}

	public bool SplitQR(Matrix mtxQ)
	{
		if (numRows < numColumns)
		{
			return false;
		}
		if (!mtxQ.Init(numRows, numRows))
		{
			return false;
		}
		for (int i = 0; i <= numRows - 1; i++)
		{
			for (int j = 0; j <= numRows - 1; j++)
			{
				int l = i * numRows + j;
				mtxQ.elements[l] = 0.0;
				if (i == j)
				{
					mtxQ.elements[l] = 1.0;
				}
			}
		}
		int nn = numColumns;
		if (numRows == numColumns)
		{
			nn = numRows - 1;
		}
		for (int k = 0; k <= nn - 1; k++)
		{
			double u = 0.0;
			int l = k * numColumns + k;
			for (int i = k; i <= numRows - 1; i++)
			{
				double w = Math.Abs(elements[i * numColumns + k]);
				if (w > u)
				{
					u = w;
				}
			}
			double alpha = 0.0;
			for (int i = k; i <= numRows - 1; i++)
			{
				double t = elements[i * numColumns + k] / u;
				alpha += t * t;
			}
			if (elements[l] > 0.0)
			{
				u = 0.0 - u;
			}
			alpha = u * Math.Sqrt(alpha);
			if (alpha == 0.0)
			{
				return false;
			}
			u = Math.Sqrt(2.0 * alpha * (alpha - elements[l]));
			if (u + 1.0 == 1.0)
			{
				continue;
			}
			elements[l] = (elements[l] - alpha) / u;
			for (int i = k + 1; i <= numRows - 1; i++)
			{
				int p = i * numColumns + k;
				elements[p] /= u;
			}
			for (int j = 0; j <= numRows - 1; j++)
			{
				double t = 0.0;
				for (int jj = k; jj <= numRows - 1; jj++)
				{
					t += elements[jj * numColumns + k] * mtxQ.elements[jj * numRows + j];
				}
				for (int i = k; i <= numRows - 1; i++)
				{
					int p = i * numRows + j;
					mtxQ.elements[p] = mtxQ.elements[p] - 2.0 * t * elements[i * numColumns + k];
				}
			}
			for (int j = k + 1; j <= numColumns - 1; j++)
			{
				double t = 0.0;
				for (int jj = k; jj <= numRows - 1; jj++)
				{
					t += elements[jj * numColumns + k] * elements[jj * numColumns + j];
				}
				for (int i = k; i <= numRows - 1; i++)
				{
					int p = i * numColumns + j;
					elements[p] -= 2.0 * t * elements[i * numColumns + k];
				}
			}
			elements[l] = alpha;
			for (int i = k + 1; i <= numRows - 1; i++)
			{
				elements[i * numColumns + k] = 0.0;
			}
		}
		for (int i = 0; i <= numRows - 2; i++)
		{
			for (int j = i + 1; j <= numRows - 1; j++)
			{
				int p = i * numRows + j;
				int l = j * numRows + i;
				double t = mtxQ.elements[p];
				mtxQ.elements[p] = mtxQ.elements[l];
				mtxQ.elements[l] = t;
			}
		}
		return true;
	}

	public bool SplitUV(Matrix mtxU, Matrix mtxV, double eps)
	{
		double[] fg = new double[2];
		double[] cs = new double[2];
		int m = numRows;
		int n = numColumns;
		if (!mtxU.Init(m, m) || !mtxV.Init(n, n))
		{
			return false;
		}
		int ka = Math.Max(m, n) + 1;
		double[] s = new double[ka];
		double[] e = new double[ka];
		double[] w = new double[ka];
		int it = 60;
		int k = n;
		if (m - 1 < n)
		{
			k = m - 1;
		}
		int l = m;
		if (n - 2 < m)
		{
			l = n - 2;
		}
		if (l < 0)
		{
			l = 0;
		}
		int ll = k;
		if (l > k)
		{
			ll = l;
		}
		if (ll >= 1)
		{
			for (int kk = 1; kk <= ll; kk++)
			{
				double d;
				if (kk <= k)
				{
					d = 0.0;
					for (int i = kk; i <= m; i++)
					{
						int ix = (i - 1) * n + kk - 1;
						d += elements[ix] * elements[ix];
					}
					s[kk - 1] = Math.Sqrt(d);
					if (s[kk - 1] != 0.0)
					{
						int ix = (kk - 1) * n + kk - 1;
						if (elements[ix] != 0.0)
						{
							s[kk - 1] = Math.Abs(s[kk - 1]);
							if (elements[ix] < 0.0)
							{
								s[kk - 1] = 0.0 - s[kk - 1];
							}
						}
						for (int i = kk; i <= m; i++)
						{
							int iy = (i - 1) * n + kk - 1;
							elements[iy] /= s[kk - 1];
						}
						elements[ix] = 1.0 + elements[ix];
					}
					s[kk - 1] = 0.0 - s[kk - 1];
				}
				if (n >= kk + 1)
				{
					for (int j = kk + 1; j <= n; j++)
					{
						if (kk <= k && s[kk - 1] != 0.0)
						{
							d = 0.0;
							for (int i = kk; i <= m; i++)
							{
								int ix = (i - 1) * n + kk - 1;
								int iy = (i - 1) * n + j - 1;
								d += elements[ix] * elements[iy];
							}
							d = (0.0 - d) / elements[(kk - 1) * n + kk - 1];
							for (int i = kk; i <= m; i++)
							{
								int ix = (i - 1) * n + j - 1;
								int iy = (i - 1) * n + kk - 1;
								elements[ix] += d * elements[iy];
							}
						}
						e[j - 1] = elements[(kk - 1) * n + j - 1];
					}
				}
				if (kk <= k)
				{
					for (int i = kk; i <= m; i++)
					{
						int ix = (i - 1) * m + kk - 1;
						int iy = (i - 1) * n + kk - 1;
						mtxU.elements[ix] = elements[iy];
					}
				}
				if (kk > l)
				{
					continue;
				}
				d = 0.0;
				for (int i = kk + 1; i <= n; i++)
				{
					d += e[i - 1] * e[i - 1];
				}
				e[kk - 1] = Math.Sqrt(d);
				if (e[kk - 1] != 0.0)
				{
					if (e[kk] != 0.0)
					{
						e[kk - 1] = Math.Abs(e[kk - 1]);
						if (e[kk] < 0.0)
						{
							e[kk - 1] = 0.0 - e[kk - 1];
						}
					}
					for (int i = kk + 1; i <= n; i++)
					{
						e[i - 1] /= e[kk - 1];
					}
					e[kk] = 1.0 + e[kk];
				}
				e[kk - 1] = 0.0 - e[kk - 1];
				if (kk + 1 <= m && e[kk - 1] != 0.0)
				{
					for (int i = kk + 1; i <= m; i++)
					{
						w[i - 1] = 0.0;
					}
					for (int j = kk + 1; j <= n; j++)
					{
						for (int i = kk + 1; i <= m; i++)
						{
							w[i - 1] += e[j - 1] * elements[(i - 1) * n + j - 1];
						}
					}
					for (int j = kk + 1; j <= n; j++)
					{
						for (int i = kk + 1; i <= m; i++)
						{
							int ix = (i - 1) * n + j - 1;
							elements[ix] -= w[i - 1] * e[j - 1] / e[kk];
						}
					}
				}
				for (int i = kk + 1; i <= n; i++)
				{
					mtxV.elements[(i - 1) * n + kk - 1] = e[i - 1];
				}
			}
		}
		int mm = n;
		if (m + 1 < n)
		{
			mm = m + 1;
		}
		if (k < n)
		{
			s[k] = elements[k * n + k];
		}
		if (m < mm)
		{
			s[mm - 1] = 0.0;
		}
		if (l + 1 < mm)
		{
			e[l] = elements[l * n + mm - 1];
		}
		e[mm - 1] = 0.0;
		int nn = m;
		if (m > n)
		{
			nn = n;
		}
		if (nn >= k + 1)
		{
			for (int j = k + 1; j <= nn; j++)
			{
				for (int i = 1; i <= m; i++)
				{
					mtxU.elements[(i - 1) * m + j - 1] = 0.0;
				}
				mtxU.elements[(j - 1) * m + j - 1] = 1.0;
			}
		}
		if (k >= 1)
		{
			for (ll = 1; ll <= k; ll++)
			{
				int kk = k - ll + 1;
				int iz = (kk - 1) * m + kk - 1;
				if (s[kk - 1] != 0.0)
				{
					if (nn >= kk + 1)
					{
						for (int j = kk + 1; j <= nn; j++)
						{
							double d = 0.0;
							for (int i = kk; i <= m; i++)
							{
								int ix = (i - 1) * m + kk - 1;
								int iy = (i - 1) * m + j - 1;
								d += mtxU.elements[ix] * mtxU.elements[iy] / mtxU.elements[iz];
							}
							d = 0.0 - d;
							for (int i = kk; i <= m; i++)
							{
								int ix = (i - 1) * m + j - 1;
								int iy = (i - 1) * m + kk - 1;
								mtxU.elements[ix] = mtxU.elements[ix] + d * mtxU.elements[iy];
							}
						}
					}
					for (int i = kk; i <= m; i++)
					{
						int ix = (i - 1) * m + kk - 1;
						mtxU.elements[ix] = 0.0 - mtxU.elements[ix];
					}
					mtxU.elements[iz] = 1.0 + mtxU.elements[iz];
					if (kk - 1 >= 1)
					{
						for (int i = 1; i <= kk - 1; i++)
						{
							mtxU.elements[(i - 1) * m + kk - 1] = 0.0;
						}
					}
				}
				else
				{
					for (int i = 1; i <= m; i++)
					{
						mtxU.elements[(i - 1) * m + kk - 1] = 0.0;
					}
					mtxU.elements[(kk - 1) * m + kk - 1] = 1.0;
				}
			}
		}
		for (ll = 1; ll <= n; ll++)
		{
			int kk = n - ll + 1;
			int iz = kk * n + kk - 1;
			if (kk <= l && e[kk - 1] != 0.0)
			{
				for (int j = kk + 1; j <= n; j++)
				{
					double d = 0.0;
					for (int i = kk + 1; i <= n; i++)
					{
						int ix = (i - 1) * n + kk - 1;
						int iy = (i - 1) * n + j - 1;
						d += mtxV.elements[ix] * mtxV.elements[iy] / mtxV.elements[iz];
					}
					d = 0.0 - d;
					for (int i = kk + 1; i <= n; i++)
					{
						int ix = (i - 1) * n + j - 1;
						int iy = (i - 1) * n + kk - 1;
						mtxV.elements[ix] = mtxV.elements[ix] + d * mtxV.elements[iy];
					}
				}
			}
			for (int i = 1; i <= n; i++)
			{
				mtxV.elements[(i - 1) * n + kk - 1] = 0.0;
			}
			mtxV.elements[iz - n] = 1.0;
		}
		for (int i = 1; i <= m; i++)
		{
			for (int j = 1; j <= n; j++)
			{
				elements[(i - 1) * n + j - 1] = 0.0;
			}
		}
		int m2 = mm;
		it = 60;
		while (true)
		{
			if (mm == 0)
			{
				ppp(elements, e, s, mtxV.elements, m, n);
				return true;
			}
			if (it == 0)
			{
				break;
			}
			int kk = mm - 1;
			while (kk != 0 && Math.Abs(e[kk - 1]) != 0.0)
			{
				double d = Math.Abs(s[kk - 1]) + Math.Abs(s[kk]);
				double dd = Math.Abs(e[kk - 1]);
				if (dd > eps * d)
				{
					kk--;
				}
				else
				{
					e[kk - 1] = 0.0;
				}
			}
			if (kk == mm - 1)
			{
				kk++;
				if (s[kk - 1] < 0.0)
				{
					s[kk - 1] = 0.0 - s[kk - 1];
					for (int i = 1; i <= n; i++)
					{
						int ix = (i - 1) * n + kk - 1;
						mtxV.elements[ix] = 0.0 - mtxV.elements[ix];
					}
				}
				for (; kk != m2 && s[kk - 1] < s[kk]; kk++)
				{
					double d = s[kk - 1];
					s[kk - 1] = s[kk];
					s[kk] = d;
					if (kk < n)
					{
						for (int i = 1; i <= n; i++)
						{
							int ix = (i - 1) * n + kk - 1;
							int iy = (i - 1) * n + kk;
							d = mtxV.elements[ix];
							mtxV.elements[ix] = mtxV.elements[iy];
							mtxV.elements[iy] = d;
						}
					}
					if (kk < m)
					{
						for (int i = 1; i <= m; i++)
						{
							int ix = (i - 1) * m + kk - 1;
							int iy = (i - 1) * m + kk;
							d = mtxU.elements[ix];
							mtxU.elements[ix] = mtxU.elements[iy];
							mtxU.elements[iy] = d;
						}
					}
				}
				it = 60;
				mm--;
				continue;
			}
			int ks = mm;
			while (ks > kk && Math.Abs(s[ks - 1]) != 0.0)
			{
				double d = 0.0;
				if (ks != mm)
				{
					d += Math.Abs(e[ks - 1]);
				}
				if (ks != kk + 1)
				{
					d += Math.Abs(e[ks - 2]);
				}
				double dd = Math.Abs(s[ks - 1]);
				if (dd > eps * d)
				{
					ks--;
				}
				else
				{
					s[ks - 1] = 0.0;
				}
			}
			if (ks == kk)
			{
				kk++;
				double d = Math.Abs(s[mm - 1]);
				double t = Math.Abs(s[mm - 2]);
				if (t > d)
				{
					d = t;
				}
				t = Math.Abs(e[mm - 2]);
				if (t > d)
				{
					d = t;
				}
				t = Math.Abs(s[kk - 1]);
				if (t > d)
				{
					d = t;
				}
				t = Math.Abs(e[kk - 1]);
				if (t > d)
				{
					d = t;
				}
				double sm = s[mm - 1] / d;
				double sm2 = s[mm - 2] / d;
				double em1 = e[mm - 2] / d;
				double sk = s[kk - 1] / d;
				double ek = e[kk - 1] / d;
				double b = ((sm2 + sm) * (sm2 - sm) + em1 * em1) / 2.0;
				double c = sm * em1;
				c *= c;
				double shh = 0.0;
				if (b != 0.0 || c != 0.0)
				{
					shh = Math.Sqrt(b * b + c);
					if (b < 0.0)
					{
						shh = 0.0 - shh;
					}
					shh = c / (b + shh);
				}
				fg[0] = (sk + sm) * (sk - sm) - shh;
				fg[1] = sk * ek;
				for (int i = kk; i <= mm - 1; i++)
				{
					sss(fg, cs);
					if (i != kk)
					{
						e[i - 2] = fg[0];
					}
					fg[0] = cs[0] * s[i - 1] + cs[1] * e[i - 1];
					e[i - 1] = cs[0] * e[i - 1] - cs[1] * s[i - 1];
					fg[1] = cs[1] * s[i];
					s[i] = cs[0] * s[i];
					if (cs[0] != 1.0 || cs[1] != 0.0)
					{
						for (int j = 1; j <= n; j++)
						{
							int ix = (j - 1) * n + i - 1;
							int iy = (j - 1) * n + i;
							d = cs[0] * mtxV.elements[ix] + cs[1] * mtxV.elements[iy];
							mtxV.elements[iy] = (0.0 - cs[1]) * mtxV.elements[ix] + cs[0] * mtxV.elements[iy];
							mtxV.elements[ix] = d;
						}
					}
					sss(fg, cs);
					s[i - 1] = fg[0];
					fg[0] = cs[0] * e[i - 1] + cs[1] * s[i];
					s[i] = (0.0 - cs[1]) * e[i - 1] + cs[0] * s[i];
					fg[1] = cs[1] * e[i];
					e[i] = cs[0] * e[i];
					if (i < m && (cs[0] != 1.0 || cs[1] != 0.0))
					{
						for (int j = 1; j <= m; j++)
						{
							int ix = (j - 1) * m + i - 1;
							int iy = (j - 1) * m + i;
							d = cs[0] * mtxU.elements[ix] + cs[1] * mtxU.elements[iy];
							mtxU.elements[iy] = (0.0 - cs[1]) * mtxU.elements[ix] + cs[0] * mtxU.elements[iy];
							mtxU.elements[ix] = d;
						}
					}
				}
				e[mm - 2] = fg[0];
				it--;
				continue;
			}
			if (ks == mm)
			{
				kk++;
				fg[1] = e[mm - 2];
				e[mm - 2] = 0.0;
				for (ll = kk; ll <= mm - 1; ll++)
				{
					int i = mm + kk - ll - 1;
					fg[0] = s[i - 1];
					sss(fg, cs);
					s[i - 1] = fg[0];
					if (i != kk)
					{
						fg[1] = (0.0 - cs[1]) * e[i - 2];
						e[i - 2] = cs[0] * e[i - 2];
					}
					if (cs[0] != 1.0 || cs[1] != 0.0)
					{
						for (int j = 1; j <= n; j++)
						{
							int ix = (j - 1) * n + i - 1;
							int iy = (j - 1) * n + mm - 1;
							double d = cs[0] * mtxV.elements[ix] + cs[1] * mtxV.elements[iy];
							mtxV.elements[iy] = (0.0 - cs[1]) * mtxV.elements[ix] + cs[0] * mtxV.elements[iy];
							mtxV.elements[ix] = d;
						}
					}
				}
				continue;
			}
			kk = ks + 1;
			fg[1] = e[kk - 2];
			e[kk - 2] = 0.0;
			for (int i = kk; i <= mm; i++)
			{
				fg[0] = s[i - 1];
				sss(fg, cs);
				s[i - 1] = fg[0];
				fg[1] = (0.0 - cs[1]) * e[i - 1];
				e[i - 1] = cs[0] * e[i - 1];
				if (cs[0] != 1.0 || cs[1] != 0.0)
				{
					for (int j = 1; j <= m; j++)
					{
						int ix = (j - 1) * m + i - 1;
						int iy = (j - 1) * m + kk - 2;
						double d = cs[0] * mtxU.elements[ix] + cs[1] * mtxU.elements[iy];
						mtxU.elements[iy] = (0.0 - cs[1]) * mtxU.elements[ix] + cs[0] * mtxU.elements[iy];
						mtxU.elements[ix] = d;
					}
				}
			}
		}
		ppp(elements, e, s, mtxV.elements, m, n);
		return false;
	}

	private void ppp(double[] a, double[] e, double[] s, double[] v, int m, int n)
	{
		int i = ((m < n) ? m : n);
		for (int j = 1; j <= i - 1; j++)
		{
			a[(j - 1) * n + j - 1] = s[j - 1];
			a[(j - 1) * n + j] = e[j - 1];
		}
		a[(i - 1) * n + i - 1] = s[i - 1];
		if (m < n)
		{
			a[(i - 1) * n + i] = e[i - 1];
		}
		for (i = 1; i <= n - 1; i++)
		{
			for (int j = i + 1; j <= n; j++)
			{
				int p = (i - 1) * n + j - 1;
				int q = (j - 1) * n + i - 1;
				double d = v[p];
				v[p] = v[q];
				v[q] = d;
			}
		}
	}

	private void sss(double[] fg, double[] cs)
	{
		double d;
		if (Math.Abs(fg[0]) + Math.Abs(fg[1]) == 0.0)
		{
			cs[0] = 1.0;
			cs[1] = 0.0;
			d = 0.0;
		}
		else
		{
			d = Math.Sqrt(fg[0] * fg[0] + fg[1] * fg[1]);
			if (Math.Abs(fg[0]) > Math.Abs(fg[1]))
			{
				d = Math.Abs(d);
				if (fg[0] < 0.0)
				{
					d = 0.0 - d;
				}
			}
			if (Math.Abs(fg[1]) >= Math.Abs(fg[0]))
			{
				d = Math.Abs(d);
				if (fg[1] < 0.0)
				{
					d = 0.0 - d;
				}
			}
			cs[0] = fg[0] / d;
			cs[1] = fg[1] / d;
		}
		double r = 1.0;
		if (Math.Abs(fg[0]) > Math.Abs(fg[1]))
		{
			r = cs[1];
		}
		else if (cs[0] != 0.0)
		{
			r = 1.0 / cs[0];
		}
		fg[0] = d;
		fg[1] = r;
	}

	public bool InvertUV(Matrix mtxAP, Matrix mtxU, Matrix mtxV, double eps)
	{
		if (!SplitUV(mtxU, mtxV, eps))
		{
			return false;
		}
		int m = numRows;
		int n = numColumns;
		if (!mtxAP.Init(n, m))
		{
			return false;
		}
		int j = n;
		if (m < n)
		{
			j = m;
		}
		j--;
		int k;
		for (k = 0; k <= j && elements[k * n + k] != 0.0; k++)
		{
		}
		k--;
		for (int i = 0; i <= n - 1; i++)
		{
			for (j = 0; j <= m - 1; j++)
			{
				int t = i * m + j;
				mtxAP.elements[t] = 0.0;
				for (int l = 0; l <= k; l++)
				{
					int f = l * n + i;
					int p = j * m + l;
					int q = l * n + l;
					mtxAP.elements[t] = mtxAP.elements[t] + mtxV.elements[f] * mtxU.elements[p] / elements[q];
				}
			}
		}
		return true;
	}

	public bool MakeSymTri(Matrix mtxQ, Matrix mtxT, double[] dblB, double[] dblC)
	{
		if (!mtxQ.Init(numColumns, numColumns) || !mtxT.Init(numColumns, numColumns))
		{
			return false;
		}
		if (dblB == null || dblC == null)
		{
			return false;
		}
		for (int i = 0; i <= numColumns - 1; i++)
		{
			for (int j = 0; j <= numColumns - 1; j++)
			{
				int u = i * numColumns + j;
				mtxQ.elements[u] = elements[u];
			}
		}
		for (int i = numColumns - 1; i >= 1; i--)
		{
			double h = 0.0;
			if (i > 1)
			{
				for (int k = 0; k <= i - 1; k++)
				{
					int u = i * numColumns + k;
					h += mtxQ.elements[u] * mtxQ.elements[u];
				}
			}
			if (h == 0.0)
			{
				dblC[i] = 0.0;
				if (i == 1)
				{
					dblC[i] = mtxQ.elements[i * numColumns + i - 1];
				}
				dblB[i] = 0.0;
			}
			else
			{
				dblC[i] = Math.Sqrt(h);
				int u = i * numColumns + i - 1;
				if (mtxQ.elements[u] > 0.0)
				{
					dblC[i] = 0.0 - dblC[i];
				}
				h -= mtxQ.elements[u] * dblC[i];
				mtxQ.elements[u] = mtxQ.elements[u] - dblC[i];
				double f = 0.0;
				for (int j = 0; j <= i - 1; j++)
				{
					mtxQ.elements[j * numColumns + i] = mtxQ.elements[i * numColumns + j] / h;
					double g = 0.0;
					for (int k = 0; k <= j; k++)
					{
						g += mtxQ.elements[j * numColumns + k] * mtxQ.elements[i * numColumns + k];
					}
					if (j + 1 <= i - 1)
					{
						for (int k = j + 1; k <= i - 1; k++)
						{
							g += mtxQ.elements[k * numColumns + j] * mtxQ.elements[i * numColumns + k];
						}
					}
					dblC[j] = g / h;
					f += g * mtxQ.elements[j * numColumns + i];
				}
				double h2 = f / (h + h);
				for (int j = 0; j <= i - 1; j++)
				{
					f = mtxQ.elements[i * numColumns + j];
					double g = (dblC[j] -= h2 * f);
					for (int k = 0; k <= j; k++)
					{
						u = j * numColumns + k;
						mtxQ.elements[u] = mtxQ.elements[u] - f * dblC[k] - g * mtxQ.elements[i * numColumns + k];
					}
				}
				dblB[i] = h;
			}
		}
		for (int i = 0; i <= numColumns - 2; i++)
		{
			dblC[i] = dblC[i + 1];
		}
		dblC[numColumns - 1] = 0.0;
		dblB[0] = 0.0;
		for (int i = 0; i <= numColumns - 1; i++)
		{
			int u;
			if (dblB[i] != 0.0 && i - 1 >= 0)
			{
				for (int j = 0; j <= i - 1; j++)
				{
					double g = 0.0;
					for (int k = 0; k <= i - 1; k++)
					{
						g += mtxQ.elements[i * numColumns + k] * mtxQ.elements[k * numColumns + j];
					}
					for (int k = 0; k <= i - 1; k++)
					{
						u = k * numColumns + j;
						mtxQ.elements[u] = mtxQ.elements[u] - g * mtxQ.elements[k * numColumns + i];
					}
				}
			}
			u = i * numColumns + i;
			dblB[i] = mtxQ.elements[u];
			mtxQ.elements[u] = 1.0;
			if (i - 1 >= 0)
			{
				for (int j = 0; j <= i - 1; j++)
				{
					mtxQ.elements[i * numColumns + j] = 0.0;
					mtxQ.elements[j * numColumns + i] = 0.0;
				}
			}
		}
		for (int i = 0; i < numColumns; i++)
		{
			for (int j = 0; j < numColumns; j++)
			{
				mtxT.SetElement(i, j, 0.0);
				switch (i - j)
				{
				case 0:
					mtxT.SetElement(i, j, dblB[j]);
					break;
				case 1:
					mtxT.SetElement(i, j, dblC[j]);
					break;
				case -1:
					mtxT.SetElement(i, j, dblC[i]);
					break;
				}
			}
		}
		return true;
	}

	public bool ComputeEvSymTri(double[] dblB, double[] dblC, Matrix mtxQ, int nMaxIt, double eps)
	{
		int m = mtxQ.GetNumColumns();
		dblC[m - 1] = 0.0;
		double d = 0.0;
		double f = 0.0;
		for (int j = 0; j <= m - 1; j++)
		{
			int it = 0;
			double h = eps * (Math.Abs(dblB[j]) + Math.Abs(dblC[j]));
			if (h > d)
			{
				d = h;
			}
			int l;
			for (l = j; l <= m - 1 && Math.Abs(dblC[l]) > d; l++)
			{
			}
			if (l != j)
			{
				do
				{
					if (it == nMaxIt)
					{
						return false;
					}
					it++;
					double g = dblB[j];
					double p = (dblB[j + 1] - g) / (2.0 * dblC[j]);
					double r = Math.Sqrt(p * p + 1.0);
					if (p >= 0.0)
					{
						dblB[j] = dblC[j] / (p + r);
					}
					else
					{
						dblB[j] = dblC[j] / (p - r);
					}
					h = g - dblB[j];
					for (int i = j + 1; i <= m - 1; i++)
					{
						dblB[i] -= h;
					}
					f += h;
					p = dblB[l];
					double e = 1.0;
					double s = 0.0;
					for (int i = l - 1; i >= j; i--)
					{
						g = e * dblC[i];
						h = e * p;
						if (Math.Abs(p) >= Math.Abs(dblC[i]))
						{
							e = dblC[i] / p;
							r = Math.Sqrt(e * e + 1.0);
							dblC[i + 1] = s * p * r;
							s = e / r;
							e = 1.0 / r;
						}
						else
						{
							e = p / dblC[i];
							r = Math.Sqrt(e * e + 1.0);
							dblC[i + 1] = s * dblC[i] * r;
							s = 1.0 / r;
							e /= r;
						}
						p = e * dblB[i] - s * g;
						dblB[i + 1] = h + s * (e * g + s * dblB[i]);
						for (int k = 0; k <= m - 1; k++)
						{
							int u = k * m + i + 1;
							int v = u - 1;
							h = mtxQ.elements[u];
							mtxQ.elements[u] = s * mtxQ.elements[v] + e * h;
							mtxQ.elements[v] = e * mtxQ.elements[v] - s * h;
						}
					}
					dblC[j] = s * p;
					dblB[j] = e * p;
				}
				while (Math.Abs(dblC[j]) > d);
			}
			dblB[j] += f;
		}
		for (int i = 0; i <= m - 1; i++)
		{
			int k = i;
			double p = dblB[i];
			if (i + 1 <= m - 1)
			{
				for (int j = i + 1; j <= m - 1 && dblB[j] <= p; j++)
				{
					k = j;
					p = dblB[j];
				}
			}
			if (k != i)
			{
				dblB[k] = dblB[i];
				dblB[i] = p;
				for (int j = 0; j <= m - 1; j++)
				{
					int u = j * m + i;
					int v = j * m + k;
					p = mtxQ.elements[u];
					mtxQ.elements[u] = mtxQ.elements[v];
					mtxQ.elements[v] = p;
				}
			}
		}
		return true;
	}

	public void MakeHberg()
	{
		int i = 0;
		for (int k = 1; k <= numColumns - 2; k++)
		{
			double d = 0.0;
			for (int j = k; j <= numColumns - 1; j++)
			{
				int u = j * numColumns + k - 1;
				double t = elements[u];
				if (Math.Abs(t) > Math.Abs(d))
				{
					d = t;
					i = j;
				}
			}
			if (d == 0.0)
			{
				continue;
			}
			if (i != k)
			{
				for (int j = k - 1; j <= numColumns - 1; j++)
				{
					int u = i * numColumns + j;
					int v = k * numColumns + j;
					double t = elements[u];
					elements[u] = elements[v];
					elements[v] = t;
				}
				for (int j = 0; j <= numColumns - 1; j++)
				{
					int u = j * numColumns + i;
					int v = j * numColumns + k;
					double t = elements[u];
					elements[u] = elements[v];
					elements[v] = t;
				}
			}
			for (i = k + 1; i <= numColumns - 1; i++)
			{
				int u = i * numColumns + k - 1;
				double t = elements[u] / d;
				elements[u] = 0.0;
				for (int j = k; j <= numColumns - 1; j++)
				{
					int v = i * numColumns + j;
					elements[v] -= t * elements[k * numColumns + j];
				}
				for (int j = 0; j <= numColumns - 1; j++)
				{
					int v = j * numColumns + k;
					elements[v] += t * elements[j * numColumns + i];
				}
			}
		}
	}

	public bool ComputeEvHBerg(double[] dblU, double[] dblV, int nMaxIt, double eps)
	{
		int n = numColumns;
		int it = 0;
		int m = n;
		while (m != 0)
		{
			int l = m - 1;
			while (l > 0 && Math.Abs(elements[l * n + l - 1]) > eps * (Math.Abs(elements[(l - 1) * n + l - 1]) + Math.Abs(elements[l * n + l])))
			{
				l--;
			}
			int ii = (m - 1) * n + m - 1;
			int jj = (m - 1) * n + m - 2;
			int kk = (m - 2) * n + m - 1;
			int ll = (m - 2) * n + m - 2;
			if (l == m - 1)
			{
				dblU[m - 1] = elements[(m - 1) * n + m - 1];
				dblV[m - 1] = 0.0;
				m--;
				it = 0;
				continue;
			}
			if (l == m - 2)
			{
				double b = 0.0 - (elements[ii] + elements[ll]);
				double c = elements[ii] * elements[ll] - elements[jj] * elements[kk];
				double w = b * b - 4.0 * c;
				double y = Math.Sqrt(Math.Abs(w));
				if (w > 0.0)
				{
					double xy = 1.0;
					if (b < 0.0)
					{
						xy = -1.0;
					}
					dblU[m - 1] = (0.0 - b - xy * y) / 2.0;
					dblU[m - 2] = c / dblU[m - 1];
					dblV[m - 1] = 0.0;
					dblV[m - 2] = 0.0;
				}
				else
				{
					dblU[m - 1] = (0.0 - b) / 2.0;
					dblU[m - 2] = dblU[m - 1];
					dblV[m - 1] = y / 2.0;
					dblV[m - 2] = 0.0 - dblV[m - 1];
				}
				m -= 2;
				it = 0;
				continue;
			}
			if (it >= nMaxIt)
			{
				return false;
			}
			it++;
			for (int j = l + 2; j <= m - 1; j++)
			{
				elements[j * n + j - 2] = 0.0;
			}
			for (int j = l + 3; j <= m - 1; j++)
			{
				elements[j * n + j - 3] = 0.0;
			}
			for (int k = l; k <= m - 2; k++)
			{
				double x;
				double y;
				double p;
				double q;
				double r;
				if (k != l)
				{
					p = elements[k * n + k - 1];
					q = elements[(k + 1) * n + k - 1];
					r = 0.0;
					if (k != m - 2)
					{
						r = elements[(k + 2) * n + k - 1];
					}
				}
				else
				{
					x = elements[ii] + elements[ll];
					y = elements[ll] * elements[ii] - elements[kk] * elements[jj];
					ii = l * n + l;
					jj = l * n + l + 1;
					kk = (l + 1) * n + l;
					ll = (l + 1) * n + l + 1;
					p = elements[ii] * (elements[ii] - x) + elements[jj] * elements[kk] + y;
					q = elements[kk] * (elements[ii] + elements[ll] - x);
					r = elements[kk] * elements[(l + 2) * n + l + 1];
				}
				if (Math.Abs(p) + Math.Abs(q) + Math.Abs(r) == 0.0)
				{
					continue;
				}
				double xy = 1.0;
				if (p < 0.0)
				{
					xy = -1.0;
				}
				double s = xy * Math.Sqrt(p * p + q * q + r * r);
				if (k != l)
				{
					elements[k * n + k - 1] = 0.0 - s;
				}
				double e = (0.0 - q) / s;
				double f = (0.0 - r) / s;
				x = (0.0 - p) / s;
				y = 0.0 - x - f * r / (p + s);
				double g = e * r / (p + s);
				double z = 0.0 - x - e * q / (p + s);
				int j;
				for (j = k; j <= m - 1; j++)
				{
					ii = k * n + j;
					jj = (k + 1) * n + j;
					p = x * elements[ii] + e * elements[jj];
					q = e * elements[ii] + y * elements[jj];
					r = f * elements[ii] + g * elements[jj];
					if (k != m - 2)
					{
						kk = (k + 2) * n + j;
						p += f * elements[kk];
						q += g * elements[kk];
						r += z * elements[kk];
						elements[kk] = r;
					}
					elements[jj] = q;
					elements[ii] = p;
				}
				j = k + 3;
				if (j >= m - 1)
				{
					j = m - 1;
				}
				for (int i = l; i <= j; i++)
				{
					ii = i * n + k;
					jj = i * n + k + 1;
					p = x * elements[ii] + e * elements[jj];
					q = e * elements[ii] + y * elements[jj];
					r = f * elements[ii] + g * elements[jj];
					if (k != m - 2)
					{
						kk = i * n + k + 2;
						p += f * elements[kk];
						q += g * elements[kk];
						r += z * elements[kk];
						elements[kk] = r;
					}
					elements[jj] = q;
					elements[ii] = p;
				}
			}
		}
		return true;
	}

	public bool ComputeEvJacobi(double[] dblEigenValue, Matrix mtxEigenVector, int nMaxIt, double eps)
	{
		int p = 0;
		int q = 0;
		if (!mtxEigenVector.Init(numColumns, numColumns))
		{
			return false;
		}
		int k = 1;
		for (int i = 0; i <= numColumns - 1; i++)
		{
			mtxEigenVector.elements[i * numColumns + i] = 1.0;
			for (int j = 0; j <= numColumns - 1; j++)
			{
				if (i != j)
				{
					mtxEigenVector.elements[i * numColumns + j] = 0.0;
				}
			}
		}
		while (true)
		{
			double fm = 0.0;
			for (int i = 1; i <= numColumns - 1; i++)
			{
				for (int j = 0; j <= i - 1; j++)
				{
					double d = Math.Abs(elements[i * numColumns + j]);
					if (i != j && d > fm)
					{
						fm = d;
						p = i;
						q = j;
					}
				}
			}
			if (fm < eps)
			{
				for (int i = 0; i < numColumns; i++)
				{
					dblEigenValue[i] = GetElement(i, i);
				}
				return true;
			}
			if (k > nMaxIt)
			{
				break;
			}
			k++;
			int u = p * numColumns + q;
			int w = p * numColumns + p;
			int t = q * numColumns + p;
			int s = q * numColumns + q;
			double x = 0.0 - elements[u];
			double y = (elements[s] - elements[w]) / 2.0;
			double omega = x / Math.Sqrt(x * x + y * y);
			if (y < 0.0)
			{
				omega = 0.0 - omega;
			}
			double sn = 1.0 + Math.Sqrt(1.0 - omega * omega);
			sn = omega / Math.Sqrt(2.0 * sn);
			double cn = Math.Sqrt(1.0 - sn * sn);
			fm = elements[w];
			elements[w] = fm * cn * cn + elements[s] * sn * sn + elements[u] * omega;
			elements[s] = fm * sn * sn + elements[s] * cn * cn - elements[u] * omega;
			elements[u] = 0.0;
			elements[t] = 0.0;
			for (int j = 0; j <= numColumns - 1; j++)
			{
				if (j != p && j != q)
				{
					u = p * numColumns + j;
					w = q * numColumns + j;
					fm = elements[u];
					elements[u] = fm * cn + elements[w] * sn;
					elements[w] = (0.0 - fm) * sn + elements[w] * cn;
				}
			}
			for (int i = 0; i <= numColumns - 1; i++)
			{
				if (i != p && i != q)
				{
					u = i * numColumns + p;
					w = i * numColumns + q;
					fm = elements[u];
					elements[u] = fm * cn + elements[w] * sn;
					elements[w] = (0.0 - fm) * sn + elements[w] * cn;
				}
			}
			for (int i = 0; i <= numColumns - 1; i++)
			{
				u = i * numColumns + p;
				w = i * numColumns + q;
				fm = mtxEigenVector.elements[u];
				mtxEigenVector.elements[u] = fm * cn + mtxEigenVector.elements[w] * sn;
				mtxEigenVector.elements[w] = (0.0 - fm) * sn + mtxEigenVector.elements[w] * cn;
			}
		}
		return false;
	}

	public bool ComputeEvJacobi(double[] dblEigenValue, Matrix mtxEigenVector, double eps)
	{
		if (!mtxEigenVector.Init(numColumns, numColumns))
		{
			return false;
		}
		for (int i = 0; i <= numColumns - 1; i++)
		{
			mtxEigenVector.elements[i * numColumns + i] = 1.0;
			for (int j = 0; j <= numColumns - 1; j++)
			{
				if (i != j)
				{
					mtxEigenVector.elements[i * numColumns + j] = 0.0;
				}
			}
		}
		double ff = 0.0;
		for (int i = 1; i <= numColumns - 1; i++)
		{
			for (int j = 0; j <= i - 1; j++)
			{
				double d = elements[i * numColumns + j];
				ff += d * d;
			}
		}
		ff = Math.Sqrt(2.0 * ff);
		ff /= 1.0 * (double)numColumns;
		bool nextLoop = false;
		while (true)
		{
			for (int i = 1; i <= numColumns - 1; i++)
			{
				for (int j = 0; j <= i - 1; j++)
				{
					double d = Math.Abs(elements[i * numColumns + j]);
					if (!(d > ff))
					{
						continue;
					}
					int p = i;
					int q = j;
					int u = p * numColumns + q;
					int w = p * numColumns + p;
					int t = q * numColumns + p;
					int s = q * numColumns + q;
					double x = 0.0 - elements[u];
					double y = (elements[s] - elements[w]) / 2.0;
					double omega = x / Math.Sqrt(x * x + y * y);
					if (y < 0.0)
					{
						omega = 0.0 - omega;
					}
					double sn = 1.0 + Math.Sqrt(1.0 - omega * omega);
					sn = omega / Math.Sqrt(2.0 * sn);
					double cn = Math.Sqrt(1.0 - sn * sn);
					double fm = elements[w];
					elements[w] = fm * cn * cn + elements[s] * sn * sn + elements[u] * omega;
					elements[s] = fm * sn * sn + elements[s] * cn * cn - elements[u] * omega;
					elements[u] = 0.0;
					elements[t] = 0.0;
					for (j = 0; j <= numColumns - 1; j++)
					{
						if (j != p && j != q)
						{
							u = p * numColumns + j;
							w = q * numColumns + j;
							fm = elements[u];
							elements[u] = fm * cn + elements[w] * sn;
							elements[w] = (0.0 - fm) * sn + elements[w] * cn;
						}
					}
					for (i = 0; i <= numColumns - 1; i++)
					{
						if (i != p && i != q)
						{
							u = i * numColumns + p;
							w = i * numColumns + q;
							fm = elements[u];
							elements[u] = fm * cn + elements[w] * sn;
							elements[w] = (0.0 - fm) * sn + elements[w] * cn;
						}
					}
					for (i = 0; i <= numColumns - 1; i++)
					{
						u = i * numColumns + p;
						w = i * numColumns + q;
						fm = mtxEigenVector.elements[u];
						mtxEigenVector.elements[u] = fm * cn + mtxEigenVector.elements[w] * sn;
						mtxEigenVector.elements[w] = (0.0 - fm) * sn + mtxEigenVector.elements[w] * cn;
					}
					nextLoop = true;
					break;
				}
				if (nextLoop)
				{
					break;
				}
			}
			if (nextLoop)
			{
				nextLoop = false;
				continue;
			}
			nextLoop = false;
			if (ff < eps)
			{
				break;
			}
			ff /= 1.0 * (double)numColumns;
		}
		for (int i = 0; i < numColumns; i++)
		{
			dblEigenValue[i] = GetElement(i, i);
		}
		return true;
	}
}
